Exact stochastic simulation of dissipation and non-Markovian effects in open quantum 

systems 



00 

o 
o 

(N 



Denis Lacroix. 

GANIL, CEA et IN2P3, BP 5027, 14076 Caen Cedex, France 
(Dated: March 31, 2008) 

The exact dynamics of a system coupled to an environment can be described by an integro- 
differential stochastic equation of its reduced density. The influence of the environment is incorpo- 
rated through a mean-field which is both stochastic and non-local in time and where the standard 
two-times correlation functions of the environment appear naturally. Since no approximation is 
made, the presented theory incorporates exactly dissipative and non-Markovian effects. Applica- 
tions to the spin-boson model coupled to a heat-bath with various coupling regimes and temperature 
show that the presented stochastic theory can be a valuable tool to simulate exactly the dynam- 
ics of open quantum systems. Links with stochastic Schroedinger equation method and possible 
extensions to "imaginary time" propagation are discussed. 
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I. INTRODUCTION 

Numerous concepts in our understanding of quantum 
mechanics have emerged from the understanding and de- 
scription of a system coupled to an environment: mea- 
surement, decoherence, appearance of classical world, ir- 
reversible process and dissipation... All these phenom- 
ena which are often encompassed in the "theory of open 
quantum systems" , bridge different fields of physics and 
chemistry [l|, 0) During the past decade, important 
advances have been made in the approximate and ex- 
act description of system embedded in an environment 
using stochastic methods. Recently the description of 
open quantum systems by Stochastic Schroedinger Equa- 
tion (SSE) has received much attention 0, 0, @- Nowa- 
days, Monte-Carlo wave-function techniques are exten- 
sively used to treat master equations in the weak coupling 
and /or Markovian limit 0, i, H, 0, II i, [H d ■ 

Large theoretical efforts are actually devoted to the 
introduction of non-Markovian effects. Among the most 
recent approaches, one could mention either determin- 
istic ap proaches like Projection Operator techniques 
0, Il4l | or stochastic methods like Quantum State Dif- 
fusion (QSD) 0, [H, [13, [11 where non-Markovian ef- 
fects are accounted for through non-local memory kernels 
and state vectors evolve according to integro-diffcrcntial 
stochastic equations. In some cases, these methods are 
shown to be exact d, [l^ . 

Recently, alternative exact methods [10, [Ml have been 
developed to treat the system-|-environment problem 
that avoid evaluation of non-local memory kernel al- 
though non-Markovian effects are accounted for exactly. 
However, up to now, only few applications of these ex- 
act techniques exist |23, [lH, [H, [2^, 2^, 25]- In all cases, 
accurate description of the short time dynamics is ob- 
tained but long time evolution can hardly be described 
due to the large increase of statistical errors with time. 
Although, the possibility to simulate exactly the dissi- 
pative dynamics of open quantum systems is already an 



important step, the challenge to describe long-time evo- 
lution is highly desirable to make the techniques more 
powerful and versatile. 



II. EXACT STOCHASTIC EQUATION FOR 
THE REDUCED SYSTEM DENSITY 

A. Introduction 

In the present work, starting from the exact stochas- 
tic formulation of ref. [2l| and projecting out the effect 
of the environment, an equation of motion for the re- 
duced system dynamics is obtained where the environ- 
ment effect is incorporated through a stochastic mean- 
field which turns out to be non-local in time. Advan- 
tages of the new stochastic theory for the description of 
long-time evolution of open quantum systems are under- 
lined. We consider here a system (S) -I- environment (E) 
described by a Hamiltonian 



H = hs + hs + hj, 



(1) 



where hs and hs denote the system and environment 
Hamiltonians respectively while hj is responsible for the 
coupling. Here we assume that the interaction writes 



ft,/ = Q(8)B, 



(2) 



where Q = /({(9Ji=i,ns) and B = g{{B^}i=i,nB) cor- 
respond to functions of two sets of operators of the sys- 
tem and environment, respectively. In particular, this 
definition includes non-linear couplings. For the sake 
of simplicity, we assume an initial separable density 
D{to) = psito) ® Psito)- As will be discussed below this 
assumption could be relaxed. The exact evolution of the 
system is described by the Liouville von-Ncumann equa- 
tion ihD — [H, D] . Due to the coupling, the simple sep- 
arable structure of the initial condition is not preserved 
in time. It has however been realized recently in sev- 
eral works using either SSE or path integral techniques 
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that the exact density of the total system D{t) could be 
obtained as an average over simple separable densities, 
i.e. D{t) — ps{t) ® PB{t)- In its simplest version, the 
stochastic process takes the form [2l| 

r dps = §-[hs,Ps] + d(,sClps + dXspsQ 
{ (3) 
I '^PE = ^[hE^PE] + d^E'BpE + dXsPs'B 

where the Ito convention for stochastic calculations is 
used [2^. d^s/E and dXg/E denote Markovian Gaussian 
stochastic variables with zero means and variances 



d^sd£.E = -T, dXsdXE = -t^, (4) 
in in 

disdXE = dXsd^E = 0. (5) 

The average over stochastic paths described by Eqs. ([3]) 
match the exact evolution. Indeed, assuming that at time 
t the total density writes D{t) = ps{t) ^ psit), the aver- 
age evolution over a small time step dt is given by 

dD = dps 'S) Pe + Ps ® dpE + dps <Si dps- (6) 

Using statistical properties of stochastic variables (Eqs. 
dUlS])), we obtain 



dps® PE + Ps®dpE = -r[hs + hE,ps® Pe], (7) 



dps®dpE = -T-[Q®^,Ps® Pe]- (8) 



Therefore, the last term simulates the interaction Hamil- 
tonian exactly and the average evolution of the total den- 
sity over a time step dt reads 

dD^^[H,D], (9) 
in 

which is nothing but the exact evolution. Here, the ex- 
actness of the method is proved assuming that the den- 
sity D{t) is a single separable density. In practice, the 
total density at time t is already an average over sepa- 
rable de nsities obtaine d along each stochastic path, i.e. 
D{t) — ps{t) ig) psit). Since Eq. ([9]) is valid for any 
density written as ps{t) ® pE{t), by summing individual 
contributions, we deduce that the evolution of the to- 
tal density obtained by averaging over different paths is 
given by ifidD = dt[H, D{ty\ which is valid at any time 
and corresponds to the exact system+cnvironment dy- 
namics. 



B. Stochastic mean- field dynamics 

Here, a slightly modified version of the stochas- 
tic process is used. It incorporates part of the sys- 
tem/environment coupling using a "mean- field" approx- 
imation in the deterministic evolution. Following ref. 
[2l| . we consider the coupled stochastic evolutions, called 
hereafter "Stochastic Mean-Field" (SMF) : 



dps = ^[hs + m))ECi,Ps]+dUQi-m))s)ps+dXsPs{Q-m))s) 



dps 



where 



§[hE + (Q(t)>sB, pe] + d^E{B ~ {B{t))E)pE + dXEPE{B - {B{t))E) 



{Qit))s = MQpsit)), {Bit))E = TT{BpE{t)). 



(10) 



(11) 



The SMF version also provides an exact reformulation of the system-j-environment problem. Indeed, the two terms 
in Eqs. ([TllHl) now read 



dps ® Pe + Ps ® dpE 



dps ® dpE 

and properly recombine to recover equation 



-T-[hs - 

in 

dt 
in 



(B(t))£;Q, PS ®pe\ + % [tlE + (Q{t))sB, PS ® pe] 
in 

{Q{t))s)®{B- {B{t))E),Ps®PE]■ 



r 



Properties of the SFM theory 



Eqs. pOj) have several advantages compared to the 
simple version (Eqs. ([3])). First, this stochastic process 
automatically insures Tt(pe) — Tt{ps) = 1 along the 
stochastic path. In addition, the inclusion of a "mean- 



field" term in the deterministic part will always reduce 
the statistical dispersion compared to the simple stochas- 
tic process given by Eqs. This reduction could be 
significant if quantum fluctuations of the coupling oper- 
ators Q and B remain small along each path 2l|, |27| . 
Indeed, at any time, a measure of the statistical fluctua- 
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tions is given by 

Xstat = Tr{ (i?t(t) _ mj)') (Z?(t) - :D(I)) } 

= Tr|L>t7:»(f)| - Tr|l3(i)^|. (12) 

Starting from the total density associated to a pure state, 
the evolution of Xstat over a small time step is directly 
connected to the average quantum fluctuations of Q and 
B, i.e. 



where, we have assumed implicitly that all second mo- 
ments except those given in equations ([4][5|) cancel out. 
For comparison, the growth of statistical fluctuations 
associated to the stochastic process without mean-field 
(Eqs. JS])) reads 

dXstat = ^{WTs+Wh}- (14) 

Eq. illustrates that the number of trajectories re- 

quired to simulate the system dynamics will depend on 



the importance of quantum fluctuations of coupling oper- 
ators along each path. In addition, a comparison between 
Eqs. (fT3|l and (fT4|) illustrates that the introduction of 
mean-field will always improve the numerical accuracy. 



C. Reduced system evolution 

In Eq. Ho]), the influence of the environment on the 
system only enters through (B(t)) e- One expects in gen- 
eral to simplify the treatment by directly considering this 
observable evolution instead of the complete pE evolu- 
tion. To express (B{t))E, we introduce the environment 
evolution operator 

UE{t,t') = exp(J^hE{t~t')^ . (15) 

Defining the new set of stochastic variables dvE and duE 
through d£,E = dvE — iduE and dXE = dvE -\- iduE, the 
evolution of pE{t) can be integrated as 



/"* ds 

PE{t) = UE{t,to)PE{to)Ul{tM)+ UE{t,s)[B,pE{s)]Ul{t,s) 

+ / dvE{s) UE{t, s){'B — {'B{s))e, Pe is)}+ulit,s)-t duEis) UEit,s)[B^{B{s)) E, Pe 
Jo Jo 

Introducing also the new variables dus and dvs defined as d^s = dus — idvs and dXs = dus + idvs, the stochastic 
equation on the reduced density reads 

dps = ^ [Hsit). ps] + dusiQ - {Qit))s,ps}+ ~ ^dvs[Q - {Cl{t))s, Ps] (17) 
in 

with Hsit) = hs + (B(i))^ Q and where the source term (B(i))^ takes the exact form 

(B(i))s = Tr{B'{t-to)pEito))-^ I D{t,s){Q{s))gds~ f D{t, s)duEis) + f D^it, s)dvE{s). (18) 

Jo Jo Jo 

I 



Here, B' {t - s) = ul{t, s)BUE{t, s) while D and Di are 
defined by: 

D{t,s) = i{[B,B'{t-s)])E, (19) 
Diit,s) = ({B-(Bis))E,B'it-s)}+)E, (20) 

where the environment expectation values are taken at 
time s, i.e. (• • •)s = Tr(- • • /9£;(s)). The two coupled 
equations, Eqs. (|17m8p provide an exact reformulation 
of the system evolution if dus/E and dvg/E verify 

dusduE = dvsdvE = ^, 
dusdvE = dvsduE = 0. 



In the following, we simply assume that the first term 
in Eq. (fT8|) cancels out. Substituting Eq. (flS)) into Eq. 
(fT7|l . we finally obtain an integro-differential stochastic 
equation for the exact system evolution where the en- 
vironment effect has been incorporated through the two 
memory functions (|19ll20p . The interesting aspect of re- 
placing Eq. (fTU|) by Eqs. (|17m8p is that, in many physical 
situations, one can generally take advantage of specific 
commutation/anti-commutation properties of B as well 
as flexibility in the noise to obtain an explicit form of the 
memory functions. 
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III. APPLICATION TO SYSTEM COUPLED TO 
A HEAT-BATH 

The method is illustrated for systems coupled to an 
environment of harmonic oscillators initially at thermal 
equilibrium and it shows that the present stochastic the- 
ory can be a valuable tool to simulate exactly open quan- 
tum systems. We take 

/.. = e(2|: + ^-"-"-") (22) 

and B = —J2n '^riXn 3- The statistical properties of 
stochastic variables dus/E and dvs/E specified above do 
not uniquely define the Wiener process. A simple pre- 
scription is to further assume 

f dusdus — dusduE — dvsdvs = dvEdvE — 0, 
\ dusdvs = dusdvE = 0. 

There are several advantages to this choice. First, 
stochastic calculus are greatly simplified. For instance, 
using standard techniques for system coupled to heat- 
bath [H, [2^ and Ito stochastic rules, shows that D and 
Di depend only on the time difference r = {t — s) and 
identify with the standard correlation functions (see 
appendix |^ : 

r+00 

D{t) ^ 2h dujJ{uj)sm{ujT), (24) 
Jo 

Di{t) ^ 2h dujJ{uj) coth{huj/2kBT) cos{ujt)(,25) 
Jo 

where 

denotes the spectral density. No approximation are made 
to obtain above equations, therefore the average over dif- 
ferent stochastic paths match the exact evolution of the 
system, including all non-Markovian effects. 

A. Equivalent Stochastic Schroedinger Equation 
formulation 

Several works, based on the influence functional 
method 0, [l^, Hol have led to similar stochastic equa- 
tions for the reduced density. For instance, authors of 
ref. [1^, [lOl use an evolution of (B(t))^ where the sec- 
ond term in Eq. ^TE\\ is absent. As demonstrated below, 
this term is of crucial importance for applications. In 
ref. 01 and in refs. [U |23, [ll] a stochastic formulation 
of the exact system-|-environment is given in terms of 
the Stochastic Schroedinger Equation (SSE) technique. 
Thanks to the additional stochastic rules Eq. (fT7)) 
has also its SSE counterpart, where system densities are 



replaced by ps = (021 and wave functions evolve ac- 
cording to 

r d|</>l) = {fHs(O+de5(Q-(Q(t))5)}|0l) 

{ (27) 
U(02| - {M{-§'Hsit) + dXsiQ~{Qit))s)} 

where the bath effect is again incorporated through the 
mean-field kernel. 

B. Application to the spin-boson model 

We illustrate the proposed technique to the spin-boson 
model. This model can be regarded as one of the sim- 
plest quantum open system coupled to a heat bath [sij 
which could not be integrated exactly. In addition, it 
has been often used as a benchmark for theories of open 
quantum systems [J, i, El [H, [H, [U, HI] . The system 
and coupling Hamiltonians are respectively chosen as 

hs — TiajQCx + fiscFz, hi — Ox 

where the {(Ji\i=x,y,z are the standard Pauli matrices. 
In the spin-boson model, the numerical solution of equa- 
tion (|17p for the system density, is equivalent to solving 
three non-linear coupled equations for the {(Ji) g- Fig. [1] 




5 , 10 15 



FIG. 1: Evolution of {o'z{t)}g (assuming {crz{0)}^ = 1) as a 
function of time obtained through the average of paths sim- 
ulated with Eq. (|17|) with the Markovian process described 
by Eqs. (|21)I23|) and memory functions given by Eq. (|24I25|I . 

We assume that J{uj) = rjui-^ — - — j. Two sets of parame- 
ters are used. In both cases, Ac = 5a;o, fcsT = 2hu!o and 
he = 0. The filled and open circles correspond respectively to 
TTTj — 0.2hLiJo and Trr; = 4huio and are compared respectively 
with the solid and dashed line obtained with the same set of 
parameters in Fig. 2 of ref. [s^. Results are obtained with 
2 X 10* trajectories. 



shows examples of dynamical evolution of {<yz{t))s ob- 
tained using Eq. (jl7p and averaging over stochastic tra- 
jectories both for weak (filled circle) and strong (open 



5 



circle) coupling. Results are compared with the Hierar- 
chical approach proposed in ref. [30| . This deterministic 
approach provides an alternative a priori exact formula- 
tion of open quantum systems dynamics and was origi- 
nally motivated by numerical difficulties encountered in 
the stochastic theory proposed in ref. |19|. Such diffi- 
culties do not occur in the present simulation and much 
less stochastic trajectories seem to be needed to accu- 
rately describe the dynamical evolution. Only 2 x lO"* 
trajectories have been used to obtained Fig. [T] leading to 
statistical errors close to zero (for comparison see discus- 
sion in [2^). The computer time for the two figures was 
less than an hour for the weak coupling case up to several 
hours for the strong coupling case on a standard personal 
computer. The difference in computing time comes from 
the fact that smaller numerical time step should be used 
as the coupling strength increases to achieve good nu- 
merical accuracy, the main difficulty being to properly 
evaluate time integrals in Eq. p^ . Denoting the time 
step by M, Mloq = 1.2 x IQ-^ and Mujq = 2.2 x 10^'' 
have been used for weak and strong coupling respectively. 

In the weak coupling case, results of our stochastic 
scheme displayed in Fig. [T] (filled circles) perfectly match 
the result of ref. [30] (sohd Hne). In contrary to ref. 
[2^ . statistical errors remain small even for long time 
evolution. The difference in numerical accuracy can be 
assigned to the second term in Eq. (jlSp which turns out 
to be crucial for numerical implementation. Stochastic 
simulations for strong coupling parameters (open circles) 
slightly differ from the results obtained with the hier- 
archical approach in ref. .'!() . The numerical conver- 
gence of the stochastic simulation presented in Fig. [T] 
has been checked. Therefore, the difference might be due 
to the fact that the numerical accuracy depends on the 
truncation scheme used in the hierarchy, even though the 
method of ref. [s^l is exact. 



C. Comparison with the Time-convolutionless 
method (TCL) 

The possibility to simulate exactly the system dynam- 
ics can also serve as a benchmark for other methods. For 
instance, we compared the exact stochastic scheme with 
the approximate Time-Convolutionless (TCL) projection 
operator method of ref. 0, [l^ . Figure [5] presents the re- 
sults of the exact stochastic simulation compared with 
the TCL2 method applied to the spin-boson model in 
ref. [13]. In this figure different cases corresponding to 
either low or high temperature regime and weak or strong 
coupling are presented. We see that the best agreement 
is obtained in the weak coupling and high temperature 
case. In general, the TCL2 method compares well with 
the exact simulation if the coupling is rather small. As 
the coupling increases (lower panels of figure [2]) , the dif- 
ference between the TCL technique and the exact method 
increases. Note that the TCL method seems to systemat- 
ically underestimate the damping. Note finally that the 



use of TCL4 equations instead of TCL2 does not improve 
the comparison. 
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FIG. 2: Average evolution of {o'x{t)}g (filled circles) and 
{ay{t))g (filled squares) as a function of time. The initial 
condition correspond to {(7x{0))g ~ 1. In all cases, hcjo ~ 
and Ac = lOe. /csT — 4he and ksT — 20he are used respec- 
tively for calculations presented in the left and right column. 
In both cases, upper panels present results for weak coupling 
{■Krj = 0.2he) while in lower panels a stronger coupling is con- 
sidered {nri — 1.0?ie). Simulations have been performed with 
4 X 10^ trajectories. For all cases, dynamical evolutions of the 
X and y spin components obtained with the TCL2 method 
0, [3 a^re displayed by solid and dashed lines respectively. 



D. Discussion of approximate system evolution 
obtained with real noise and Hermitian system 
densities in stochastic evolution 

Applications presented above use specific constraints 
on the Markovian process given by Eqs. This pre- 

scription greatly simplifies stochastic calculus. For in- 
stance, simple exact expressions have been obtained for 
D{t, s) and Di{t, s) when the system is coupled to a heat- 
bath of harmonic oscillators (Eqs. (|24II25[) ). The main 
consequence of Eqs. (j23p is that dug/E ^-nd dLvg/^ should 
be complex stochastic variables leading to non Hermi- 
tian densities along paths. As illustrated above, such 
a stochastic process could be used to simulate exactly 
the reduced density evolution. The main disadvantage 
of non-Hermitian densities is however that system ob- 
servables could hardly be interpreted. We discuss here 
the possibility to perform stochastic evolution of reduced 
Hermitian densities. 

Relaxing the constraints given by Eqs. (|23p . autho- 
rizes to choose dus/E and dvs/E as real stochastic vari- 
ables, which automatically insures that ps{t) and pE(t) 
remain hermitian. This alternative has however two ma- 
jor drawbacks. First, one cannot anymore have an equiv- 
alent SSE picture. Second, while D{t, s) still identi- 



6 



fies with Eq. p4)) . no simple exact expression can be 
worked out for Di{t,s). However, since this kernel is a 
functional of pE^s), a hierarchy of more and more ac- 
curate approximations could be obtained by successive 
replacements of pe{s) into Eq. p9)) by its integral ex- 
pression, Eq. ()16|) . In the present work, we concen- 
trate on the simplest case where pe{s) is replaced by 
Pe{s) — UE{t,to)pE{to)Ul,{t,to) in the time integral of 
the memory kernel. In this limit, Di{s, t) also reduces to 
Eq. By doing this approximation, the stochastic 

process is not exact anymore. Fig. [3] presents a com- 
parison of the exact stochastic simulation obtain with 
complex noise (filled circles) and the approximate case 
with real noise (open circles). The parameters of the 
spin-boson model are the same as in Fig. [1] This fig- 
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FIG. 3: Results obtained in the weak and strong coupling case 
with the exact stochastic simulation using complex noises and 
non-hermitian (filled circles) are compared with the approxi- 
mate simulation (open circles) using real noise and hermitian 
densities along each paths. In this figure, the values of pa- 
rameters are the same as in Fig. [1] 



ure shows that the approximate scheme with real noise 
is very close to the exact simulation even in the strong 
coupling limit. In the latter case, only at very large time, 
the two simulation start to deviate slightly from one an- 
other. For completeness, approximate stochastic simu- 
lations obtained for cases presented in figure [2] are com- 
pared to the exact scheme in figure [H We see that except 
for the weak coupling and low temperature case, the ap- 
proximate simulation is very close to the exact case. It is 
worth mentioning that approximate simulation presented 
here uses the simplest prescription for Di{s,t). There- 
fore, improved description could a priori be obtained us- 
ing better approximation of Di{s,t) obtained with the 
method described above. This example is very encourag- 
ing and provides a new method to simulate open systems 
with a stochastic process preserving the hermitian prop- 
erties of the system density. 



FIG. 4: Average evolution of {a-x{t)) g (filled circles) and 
{o'y{t))g (filled squares) as a function of time obtained with 
the exact stochastic scheme using the same sets of parameters 
as in figure[2] In each cases, dynamical evolutions of the x and 
y spin components obtained with the approximate stochas- 
tic simulation using real noise are respectively displayed with 
open circles and open squares. 



IV. CONCLUSION 

The results obtained with the stochastic theory for 
the spin-boson model are very encouraging. The the- 
ory turns out to be accurate not only to simulate short 
but also long-time dynamics and does not seem to suf- 
fer from numerical instability quoted in ref 2J|. It is 
worth mentioning that optimization techniques proposed 
in ref. 21] which are not used here, can be implemented 
to further reduce the number of stochastic paths. Besides 
the numerical aspects, the link with classical dissipations 
dynamics could be easily made, similarly to ref. 

Our stochastic approach could be generalized to any 
initial conditions that could be written as a mixing of 
separable densities, i.e. D{to) = X^n (*o) with 

= p's ® p'^E where the VF„ are complex coefficients. 
Then, the complete exact dynamics is recovered by both 
averaging over trajectories starting from each (t = 0) 
individually and averaging over the initial conditions. 
The theory is also not restricted to real time evolution. 
Statistical properties of the system-|-environment could 
be studied by considering imaginary time propagation, 
i.e. idt/h — > (3. Imaginary time propagations lead nat- 
urally to densities written as a mixing of separable den- 
sities and could then serve as initial conditions for real 
time evolution. By combining both imaginary time and 
real time propagation, general physical problems similar 
to those depicted in ref. @ can be addressed. The main 
limitation of the technique is clearly the choice of cou- 
pling operator B which should give simple memory func- 
tions (Eq. - 1^0)) ). It is however worth mentioning 
that most of the coupling operators used in the context 
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of open quantum systems already enter into this cate- 
gory [l|, 0, [13] ■ We believe that the stochastic theory 
presented in this paper can be a valuable tool to access 
exactly the dynamics of more complex open quantum sys- 
tems. We presented here specific applications on systems 
coupled to a heat-bath. The approach can however be 
applied to various types of environments and couplings 
which might be of great interest to address dissipation, 
measurement and/or decoherence problems in quantum 
systems. 
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APPENDIX A: PROOF OF EQ. ^ FOR A 
HEAT-BATH OF HARMONIC OSCILLATORS 

We give here a proof of Eq. (fTS]) where D and Di iden- 
tify with Eqs. (|24j|25p . The environment is assumed to 
be a set of harmonic oscillators, labeled by "n" associated 
to creation/annihilation (ajj,a„), i.e. 

H ^^huJnUian + ^j . (Al) 



d{an) 
d{al) 

where we have introduced the notation c„ = ~Kn/ \/2'i]n 
and rjn = m„a;„/?i. Here, cr^|"_^, cri"l and (7^*2 denote the 
second moments of the aj^, a„ operators: 

cr+2 = (4a„) - (a„) (at ) = ctL"! - 1, 
cri"2 = (a„a„) - (a„) (a„) , 

cr+l = (44) - (4) (4) • 

According to the stochastic environment dynamics, we 
can show that these moments simply evolve as 

' d4"2 = -2iw„di4"l, 

< d4"| = +2ic^„di4"|, (A4) 



If thermal equilibrium is initially assumed, the environ- 
ment density writes as a product of densities of each os- 
cillator, pe = TlnPn where each density can be written 
as Gaussian operators (see for instance 1^) determined 
by first and second moments of the (at,a„). The time 
evolution of the environment is given by Eq. ()10|) where 
B = — Kn^n and where the fiuctuating variables ver- 
ify (according to Eq. (^5)) ): 



(IXecIXe ^ d^Ed(,E = dXEd^E = 0. (A2) 



Above prescription and the specific form of B induces 
important simplifications listed below. First, the ini- 
tial product form of the environment density is preserved 
along the stochastic paths where each oscillator density 
verifies Tr(p„) — 1 and where for each pairs of densities 
d{pnPm) " Pndprn + dpnPm- Secoud, duc to the linear 
coupling operator, the Gaussian nature of initial den- 
sities is also preserved along paths. Therefore, we can 
equivalently solve the density equation of motion or fol- 
low first and second moments of each (at,a„) in time. 
Here we consider the second strategy. 

From the pE evolution, the equation of motion of each 
pairs (at ) = Tr(p£;at ) and (a„) = Tr(p£;a„) reads: 



Since we assume that each oscillator is initially at ther- 
mal equilibrium, we deduce that second moments are 
constant in time with a^^l{t) = a^_^^{t) — 0, while 
CT-Iit) = cr+2 + 1 = iN{ujn) + 1). Here, we have in- 
troduced the standard function [12]: (2iV(ijj„) -f 1) = 
coth(a[j„/(2fcBT)). Substituting in Eqs. (|A3iA4p and 
using standard integration techniques ^ leads finally to 



n ^ ^ 
I 

= -iuj^dt (a„) + ^cn (Q) -I- c„ {d^E (t) + (tL"2 (t) -f l) -f dXE [<J^+1 it) + a^l'l (i)) } , (A3) 

- +zc.„dt (at ) - ^c„ (Q) + c„ [d^E (^l"2(t) + + dXE (^i"2(t) + 4"j(i) + l) } , 

I 



daX^' = 0. 



(a„(0) = e-''^"* (a„(0)) + ^ f ^~^^,.(t-s) (q^^)) + j' |^^^(^) ^^^^^^ ^ ^ dXE{s)N{u^,,)] , 

JO JO 
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(at (i)) = e+-"* (at (0)) f ^+^^^(t-s) ^q(^)^ ^ T + dXE{s){N{u:,,) + 1)}(.A5) 

Accordingly, each position operator entering in (B(<)) reads: 

= -^((a„(0))e— "* + (at(0))e+^-"*| + ^ T sin(u;„[t - s]) (Q(s)) ds, 

-:T^ / {cosK[t - s])(27VK) + 1) [d^E{s) + dA£(s)] - zsinK[i - s\) [d^E^s) - dA^Gs)]} . 

-^"nn Jo 

Assuming the initial conditions (a„(0)) = (at(0)) = 0, substituting in the expression of (B(i)) — ~ X)n '^n (^^"(O) 
and introducing the spectral density (Eq. (pS)) ). we finally recover equation where £> and Di are respectively 
given by Eqs. dM]) and ((25)) . 
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